home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / dgamma.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  6.4 KB  |  131 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((ngam 0)
  12.       (xmin 0.0)
  13.       (xmax 0.0)
  14.       (dxrel 0.0)
  15.       (gamcs (make-array 42 :element-type 'double-float))
  16.       (pi_ 3.141592653589793)
  17.       (sq2pil 0.9189385332046728)
  18.       (first nil))
  19.   (declare (type f2cl-lib:logical first)
  20.            (type (simple-array double-float (42)) gamcs)
  21.            (type double-float sq2pil pi_ dxrel xmax xmin)
  22.            (type f2cl-lib:integer4 ngam))
  23.   (f2cl-lib:fset (f2cl-lib:fref gamcs (1) ((1 42))) 0.00857119559098933)
  24.   (f2cl-lib:fset (f2cl-lib:fref gamcs (2) ((1 42))) 0.004415381324841007)
  25.   (f2cl-lib:fset (f2cl-lib:fref gamcs (3) ((1 42))) 0.05685043681599364)
  26.   (f2cl-lib:fset (f2cl-lib:fref gamcs (4) ((1 42))) -0.00421983539641856)
  27.   (f2cl-lib:fset (f2cl-lib:fref gamcs (5) ((1 42))) 0.0013268081812124603)
  28.   (f2cl-lib:fset (f2cl-lib:fref gamcs (6) ((1 42))) -1.8930245297988807e-4)
  29.   (f2cl-lib:fset (f2cl-lib:fref gamcs (7) ((1 42))) 3.606925327441246e-5)
  30.   (f2cl-lib:fset (f2cl-lib:fref gamcs (8) ((1 42))) -6.056761904460864e-6)
  31.   (f2cl-lib:fset (f2cl-lib:fref gamcs (9) ((1 42))) 1.0558295463022833e-6)
  32.   (f2cl-lib:fset (f2cl-lib:fref gamcs (10) ((1 42))) -1.811967365542384e-7)
  33.   (f2cl-lib:fset (f2cl-lib:fref gamcs (11) ((1 42))) 3.117724964715322e-8)
  34.   (f2cl-lib:fset (f2cl-lib:fref gamcs (12) ((1 42))) -5.354219639019687e-9)
  35.   (f2cl-lib:fset (f2cl-lib:fref gamcs (13) ((1 42))) 9.193275519859591e-10)
  36.   (f2cl-lib:fset (f2cl-lib:fref gamcs (14) ((1 42))) -1.57794128028834e-10)
  37.   (f2cl-lib:fset (f2cl-lib:fref gamcs (15) ((1 42))) 2.7079806229349546e-11)
  38.   (f2cl-lib:fset (f2cl-lib:fref gamcs (16) ((1 42))) -4.64681865382573e-12)
  39.   (f2cl-lib:fset (f2cl-lib:fref gamcs (17) ((1 42))) 7.97335019200742e-13)
  40.   (f2cl-lib:fset (f2cl-lib:fref gamcs (18) ((1 42))) -1.368078209830916e-13)
  41.   (f2cl-lib:fset (f2cl-lib:fref gamcs (19) ((1 42))) 2.3473194865638006e-14)
  42.   (f2cl-lib:fset (f2cl-lib:fref gamcs (20) ((1 42))) -4.027432614949067e-15)
  43.   (f2cl-lib:fset (f2cl-lib:fref gamcs (21) ((1 42))) 6.910051747372101e-16)
  44.   (f2cl-lib:fset (f2cl-lib:fref gamcs (22) ((1 42))) -1.185584500221993e-16)
  45.   (f2cl-lib:fset (f2cl-lib:fref gamcs (23) ((1 42))) 2.034148542496374e-17)
  46.   (f2cl-lib:fset (f2cl-lib:fref gamcs (24) ((1 42))) -3.490054341717406e-18)
  47.   (f2cl-lib:fset (f2cl-lib:fref gamcs (25) ((1 42))) 5.987993856485305e-19)
  48.   (f2cl-lib:fset (f2cl-lib:fref gamcs (26) ((1 42))) -1.027378057872228e-19)
  49.   (f2cl-lib:fset (f2cl-lib:fref gamcs (27) ((1 42))) 1.76270281606053e-20)
  50.   (f2cl-lib:fset (f2cl-lib:fref gamcs (28) ((1 42))) -3.0243206537353057e-21)
  51.   (f2cl-lib:fset (f2cl-lib:fref gamcs (29) ((1 42))) 5.188914660218398e-22)
  52.   (f2cl-lib:fset (f2cl-lib:fref gamcs (30) ((1 42))) -8.902770842456577e-23)
  53.   (f2cl-lib:fset (f2cl-lib:fref gamcs (31) ((1 42))) 1.527474068493343e-23)
  54.   (f2cl-lib:fset (f2cl-lib:fref gamcs (32) ((1 42))) -2.620731256187363e-24)
  55.   (f2cl-lib:fset (f2cl-lib:fref gamcs (33) ((1 42))) 4.496464047830538e-25)
  56.   (f2cl-lib:fset (f2cl-lib:fref gamcs (34) ((1 42))) -7.714712731336878e-26)
  57.   (f2cl-lib:fset (f2cl-lib:fref gamcs (35) ((1 42))) 1.3236354531260444e-26)
  58.   (f2cl-lib:fset (f2cl-lib:fref gamcs (36) ((1 42))) -2.2709994129429292e-27)
  59.   (f2cl-lib:fset (f2cl-lib:fref gamcs (37) ((1 42))) 3.896418998003992e-28)
  60.   (f2cl-lib:fset (f2cl-lib:fref gamcs (38) ((1 42))) -6.685198115125953e-29)
  61.   (f2cl-lib:fset (f2cl-lib:fref gamcs (39) ((1 42))) 1.1469986631400242e-29)
  62.   (f2cl-lib:fset (f2cl-lib:fref gamcs (40) ((1 42))) -1.9679385863451343e-30)
  63.   (f2cl-lib:fset (f2cl-lib:fref gamcs (41) ((1 42))) 3.3764488165853374e-31)
  64.   (f2cl-lib:fset (f2cl-lib:fref gamcs (42) ((1 42))) -5.793070335782136e-32)
  65.   (setq first f2cl-lib:%true%)
  66.   (defun dgamma (x)
  67.     (declare (type double-float x))
  68.     (prog ((sinpiy 0.0) (y 0.0) (dgamma 0.0) (i 0) (n 0))
  69.       (declare (type f2cl-lib:integer4 n i)
  70.                (type double-float dgamma y sinpiy))
  71.       (cond
  72.        (first
  73.         (setf ngam
  74.                 (initds gamcs 42
  75.                  (* 0.1f0 (f2cl-lib:freal (f2cl-lib:d1mach 3)))))
  76.         (multiple-value-bind
  77.             (var-0 var-1)
  78.             (dgamlm xmin xmax)
  79.           (declare (ignore))
  80.           (setf xmin var-0)
  81.           (setf xmax var-1))
  82.         (setf dxrel (f2cl-lib:fsqrt (f2cl-lib:d1mach 4)))))
  83.       (setf first f2cl-lib:%false%)
  84.       (setf y (coerce (abs x) 'double-float))
  85.       (if (> y 10.0) (go label50))
  86.       (setf n (f2cl-lib:int x))
  87.       (if (< x 0.0) (setf n (f2cl-lib:int-sub n 1)))
  88.       (setf y (- x n))
  89.       (setf n (f2cl-lib:int-sub n 1))
  90.       (setf dgamma (+ 0.9375 (dcsevl (- (* 2.0 y) 1.0) gamcs ngam)))
  91.       (if (= n 0) (go end_label))
  92.       (if (> n 0) (go label30))
  93.       (setf n (f2cl-lib:int-sub n))
  94.       (if (= x 0.0) (xermsg "SLATEC" "DGAMMA" "X IS 0" 4 2))
  95.       (if (and (< x 0.0f0) (= (- (+ x n) 2) 0.0))
  96.           (xermsg "SLATEC" "DGAMMA" "X IS A NEGATIVE INTEGER" 4 2))
  97.       (if
  98.        (and (< x -0.5) (< (abs (/ (- x (f2cl-lib:aint (- x 0.5))) x)) dxrel))
  99.        (xermsg "SLATEC" "DGAMMA"
  100.         "ANSWER LT HALF PRECISION BECAUSE X TOO NEAR NEGATIVE INTEGER" 1 1))
  101.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  102.                     ((> i n) nil)
  103.         (tagbody (setf dgamma (/ dgamma (- (+ x i) 1))) label20))
  104.       (go end_label)
  105.      label30
  106.       (f2cl-lib:fdo (i 1 (f2cl-lib:int-add i 1))
  107.                     ((> i n) nil)
  108.         (tagbody (setf dgamma (* (+ y i) dgamma)) label40))
  109.       (go end_label)
  110.      label50
  111.       (if (> x xmax) (xermsg "SLATEC" "DGAMMA" "X SO BIG GAMMA OVERFLOWS" 3 2))
  112.       (setf dgamma 0.0)
  113.       (if (< x xmin)
  114.           (xermsg "SLATEC" "DGAMMA" "X SO SMALL GAMMA UNDERFLOWS" 2 1))
  115.       (if (< x xmin) (go end_label))
  116.       (setf dgamma
  117.               (exp
  118.                (+ (- (* (- y 0.5) (f2cl-lib:flog y)) y) sq2pil (d9lgmc y))))
  119.       (if (> x 0.0) (go end_label))
  120.       (if (< (abs (/ (- x (f2cl-lib:aint (- x 0.5))) x)) dxrel)
  121.           (xermsg "SLATEC" "DGAMMA"
  122.            "ANSWER LT HALF PRECISION, X TOO NEAR NEGATIVE INTEGER" 1 1))
  123.       (setf sinpiy (sin (* pi_ y)))
  124.       (if (= sinpiy 0.0)
  125.           (xermsg "SLATEC" "DGAMMA" "X IS A NEGATIVE INTEGER" 4 2))
  126.       (setf dgamma (/ (- pi_) (* y sinpiy dgamma)))
  127.       (go end_label)
  128.      end_label
  129.       (return (values dgamma nil)))))
  130.  
  131.